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RADFLO PHYSICS AND ALGORITHMS 


by 


Eugene M. D. Symbalisty, John Zinn, and Rodney W. Whitaker 


ABSTRACT 
This paper describes the history, physics, and algorithms of the 
computer code RADFLO and its extension HYCHEM. RADFLO is a 
one-dimensional, radiation-transport hydrodynamics code that is used 
to compute early-time fireball behavior for low-altitude nuclear 
bursts. The primary use of the code is the prediction of optical signals 
produced by nuclear explosions. It has also been used to predict thermal 
and hydrodynamic effects that are used for vulnerability and lethality 
applications. Another closely related code, HYCHEM, is an extension 
of RADFLO which includes the effects of nonequilibrium 
chemistry. Some examples of numerical results will be shown, along 
with scaling expressions derived from those results. We describe new 
computations of the structures and luminosities of steady-state shock 
waves and radiative thermal waves, which have been extended to cover 
a range of ambient air densities for high-altitude applications. We also 
describe recent modifications of the codes to use a one-dimensional 
analog of the CAVEAT fluid-dynamics algorithm in place of the former 
standard Richtmyer-von Neumann algorithm. 
_______________________________________ 


1.0 INTRODUCTION 
RADFLO is a spherically-symmetric, multifrequency, radiation-transport 
and Lagrangian hydrodynamics code designed for modeling nuclear fireballs in the 
atmosphere for altitudes up to about 50 km. The model assumes local thermodynamic 
equilibrium and makes use of stored tables of the air equation of state (pressure, temperature, 
gamma, and sound speed as functions of density and internal energy) and of air opacity (vs 
density, temperature, and photon energy). A related code, HYCHEM, includes RADFLO 
but adds nonequilibrium chemistry and effects of gamma rays and neutrons as well. 
As the code runs it generates profiles of temperature, density, material 
velocity, pressure, and brightness at various wavelengths—for a sequence of times. At the 
end of the run, RADFLO produces plots of radiated power vs time in several wavelength 
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bands along with plots of fireball radius vs time and expansion velocity vs time. RADFLO 
is also used to provide initial conditions for related problems (Jones et al. 1987, Symbalisty 
et al. 1987, Symbalisty 1991). 
The code is designed to treat the evolution of the air fireball as accurately as possible. 
The bomb itself and associated bomb debris are treated relatively crudely. The debris 
equation of state and opacities are taken to be the same as those of air with a few exceptions, 
which are described later in this report. Also, we do not attempt to compute the Taylor- 
unstable debris expansion; the debris cloud radius is modeled with a crude empirical 
algorithm. 
The altitude range of applicability of the code is limited mainly by the assumption 
of spherical symmetry. That is, the radiation transport and hydrodynamic motions must 
produce a spherical fireball approximately concentric with the bomb. This implies roughly 
that, for the dominant wavelengths of radiation emitted by the bomb, the mean free paths 
must be small compared to the atmospheric scale height. Also the fireball radius (during 
times of interest) must be smaller than the scale height. For a Megaton device this implies 
an altitude limit of about 50 km; however, for a small yield device with few energetic x-rays 
the limit can be higher. The term yield is the sum of the energy radiated from the case in x- 
rays, the residual debris’ internal energy, and the debris’ kinetic energy. 
A distinctive feature of RADFLO is the algorithm used to compute the radiation 
transport. Functionally it is a “two-stream” approximation, implying that for each frequency 
group the intensity function has one value for polar angles smaller than a prescribed angle 
θs and another value for θ ≥ θs. The angle θs is computed at each radial location and time 
and corresponds to the angle subtended by the opaque “fireball.” Inside the fireball itself, θs 
is set to π / 2. In regions that are optically thick, the diffusion approximation is used. The 
transport computation is done for each of 42 frequency groups ranging from the far infrared 
up through the soft x-ray region (Zinn 1973). 
A bothersome feature of the physical and numerical problem is a very bad mismatch 
of length scales between the fireball itself, which has dimensions of meters to tens of meters, 
and the thickness of the fireball surface region that determines its brightness (i.e., a shock 
front or thermal-wave front), which is typically millimeters to centimeters. To cope with 
this problem a separate computer code, SSW, computes the temperature structures and radiant 
fluxes in steady-state shock waves and thermal waves (i.e. waves that propagate with a 
time-invariant structure, such that the time can be eliminated from the radiation and 
hydrodynamic equations). The results of these steady-state computations are used 
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within RADFLO for computing the radiation fluxes emanating from the fireball front at 
early times (for wavelengths longer than 180 nm and times when the front velocity exceeds 
30 km/s). 
Over the years a number of people have worked with RADFLO at Los Alamos 
National Laboratory (LANL) or Mission Research Corporation (MRC), and possibly other 
institutions. In this report, especially in the historical sections, the present authors use the 
term “we,” which in some cases refers to these past participants. Thus “we” should be 
interpreted as the present authors and past contributors. In particular this report draws heavily 
on the description of RADFLO in four reports by Sappenfield (1973, 1974, 1975, 1976). 


2.0 RADFLO HISTORY 
The first version of RADFLO was written in 1963 shortly after the high-altitude 
nuclear test series of Operation Dominic in order to model and analyze the Bluegill event (at 
“intermediate” altitude). For that objective it was written initially as a radiation transport- 
only code without hydrodynamics (Zinn 1964). That initial version was fairly successful in 
matching the fireball optical output vs time to the end of the main optical 
pulse. Hydrodynamics was added to the code in 1964. With the hydrodynamics the Bluegill 
calculations were extended to later times (several seconds), and comparisons with the Bluegill 
measurements of brightness, radiated power, and fireball dimensions were amazingly 
good (Zinn and Fajen 1964; Zinn 1965). 
During the following year RADFLO computations were run for the Tightrope event 
(somewhat smaller yield and lower in altitude). Based on comparisons of the computed 
results with Tightrope radius-time data, the estimated Tightrope yield was revised upward. 
After an hiatus of five years due to shifts in military interest to higher altitude, 
RADFLO was again used in the early 1970s for the nuclear tests at sea level, in particular 
the French tests (Zinn et al. 1970). RADFLO did a good job of modeling a sea level explosion 
in the time regime from shortly after first maximum until shortly after second maximum, and 
was also useful for the early-time pre-first-maximum phase of large-yield explosions for 
studies of the x-ray veil (Zinn et al. 1972; Zinn, Kodis, and Sutherland 1974). However, for 
the time span between veil breakthrough and first maximum the computed brightness and 
radiated power results exhibited very large oscillations. These oscillations could be attributed 
to the necessarily coarse spatial mesh used in the code, whose dimensions were several 
orders of magnitude larger than those expected of the fireball front (i.e. a shock front or 
radiative thermal-wave front). 
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Since the detailed temperature structure of the front determines its brightness, we 
embarked on a separate study of the structure of high-velocity shock waves which included 
radiation effects. For this study we transformed the radiation and hydrodynamic equations 
into a coordinate system moving at the shock velocity, with the assumption that in the moving 
coordinate system the temperature and density structure would be time-invariant. The 
finite difference representations of these equations were incorporated into a code called 
SSS (steady-state shock). The results showed that for shock waves with velocities in the 
range 30 to 80 km/s, the brightness-determining structure had dimensions of a millimeter or 
so, and the effective brightness temperature for velocities above 30 km/s was significantly 
smaller than the shock Rankine-Hugoniot temperature. Above 80 km/s there was an apparent 
breakdown of the steady-state assumption. The SSS results were published (Zinn and 
Anderson 1973) but not incorporated into RADFLO until several years later. 
In the mid-1970s several projects compared RADFLO computations with old US 
nuclear test data and with data from the French tests (e.g. Horak et al. 1983). We studied 
relationships between x-ray veil breakthrough times and expected nuclear device x-ray outputs 
(Zinn et al. 1970), and we compared fireball outputs at X- and K-band microwave 
wavelengths with microwave radiometer data (Zinn et al. 1974). We looked at the contribution 
due to device neutrons and gamma rays of air fluorescence to the total optical power vs 
time. We made detailed comparisons of computed minimum and second maximum times 
with yield-scaling laws derived from nuclear test data (Zinn et al. 1972; Zinn, Kodis, and 
Sutherland 1974). We commenced studies on the effects of nuclear “smog” associated with 
nonequilibrium chemistry ( NO2 , ozone, etc.) outside the fireball. Also during the 1970s 
several two-dimensional late-time fireball codes were developed which used RADFLO results 


as input. 
Also in the 1970s the code began to proliferate into several separate versions. In 
1971 the then-current code was taken to MRC, where RADFLO was modified for a smaller 
computer and subsequently further modified in other ways. The name RADFLO was 
retained. The MRC RADFLO code was used extensively for analysis of satellite bhangmeter 
data and for Defense Nuclear Agency (DNA) systems studies (Sappenfield 1973, 1974, 
1975, and 1976). The MRC code was supplied under contract to the Air Force. Also, 
a reorganization within Los Alamos led to two separate RADFLO versions in two groups, 
which gradually diverged from one another. For the version upon which this report is based, 
modifications and improvements made during the 1970s were usually discussed with MRC, 
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so there are strong resemblances between this code and the MRC version. In this period we 
became interested in nonequilibrium chemistry effects and “smog,” and we continued to 
look for solutions to the first-maximum radiant-power oscillation problem. 
Only a few comments will be made on the second Los Alamos version of the 
RADFLO code mentioned in the paragraph above. This version was used for low-altitude 
(< 30 km) bursts and assumes equilibrium chemistry in the equation of state. It has been 
used for prediction of power-time curves and airblast effects, as well as initialization of 
two-dimensional calculations. Documentation can be found in Horak and Kodis (1983). 
The appendix of the Horak and Kodis report describes the use of fine zoning to achieve 


overall 3% accuracy in shock over pressures using initial zoning of 2.5cm / kt1/3. The 
report also includes a description of the use of discrete ordinate radiation-transport algorithm. 
This version has been used to calculate the evolution of concentric nuclear explosions; 
bursts at the same place but separated in time (Jones et al. 1983). 
Another area of interest has been the study of mass effects on optical-power time 
histories. The original RADFLO essentially assumed a point explosion with a very high 
yield to mass ratio. Realistic mass values were modeled with a heavy-air approximation 
(Horak 1980; Jones et al. 1980). This approximation assumes that the mass of the device or 
the mass immediately surrounding the device can be modeled, to first order, as high-density 
or heavy air. Extreme effects of shielding can be found in Symbalisty (1994) where a different 
material, with its own equation of state and opacity, is used to model the mass surrounding 
the device. The key result of large amounts of mass surrounding the device is the distortion 
of the power-time curve up to and around first maximum. This is the time regime where 
RADFLO is least equipped to give reliable answers due to the oscillation problem and the 
simple device description. 
From 1979 to 1982 we extended the steady-state shock wave model to include 
radiation transport-dominated waves (thermal waves), which led to a steady-state thermal- 
wave computer model, SST, with strong similarities to the SSS model. The two computer 
codes were merged into a single steady-state shock and thermal-wave code, SSW (Zinn and 
Sutherland 1981). The results showed that for both “shock waves” (i.e. waves dominated by 
hydrodynamics—with propagation velocities below 80 km/s) and “thermal waves” (i.e. waves 
dominated by radiation transport—with velocities above 80 km/s) the thickness of the visible 
surface layer was a centimeter or less, and the effective visible brightness temperature was 
much smaller than the actual temperature behind the wave front. The computed luminosity 
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results could be summarized in the form of plots of effective brightness temperature (for 
specified wavelength bands) vs propagation velocity, for given values of the ambient 
air density. Very recently (1995) we have extended the SSW computations to include a range 
of values of the background air density so that the results now apply to high altitudes as well 
as to sea level. These results are summarized in Appendix A of this report. 
The SSW results showed that the thickness of the fireball radiating-surface layer 
was much smaller than the smallest possible zone thicknesses that could be used in a 
functioning explicit Lagrangian fireball rad-hydro code. Therefore we needed to devise a 
way of incorporating the results of the SSW computations into the RADFLO radiation 
transport algorithm. The best approach was to compute the radiant intensity emanating from 
the fireball surface from the computed velocity of that surface, using a stored tabulation of 
brightness temperature vs velocity. The resulting RADFLO was free of most of the 
troublesome oscillations (Zinn and Sutherland 1982). However, persistent nagging questions 
remained concerning precisely where to define the fireball front, how to compute its velocity, 
and how to correct for the effects of the artificial viscous pressure. 
As noted earlier, the LANL (Zinn, Kodis, and Sutherland) and MRC (Sappenfield 
and McCartor) versions of RADFLO began to diverge in 1971. The latter version was brought 
back to LANL by D. Sappenfield in 1980 (Sappenfield 1983, 1984, and 1985), and since 
that time both codes have been maintained separately. This leads to some confusion which 
we have not quite sorted out. For now we will refer to the two versions as RADFLOz and 
RADFLOs, respectively. A vintage 1979 version of RADFLOz was used for the Vela Alert- 
747 computations and later described in a RADFLO User’s Manual (Horak 1980; Jones et 
al. 1987; Horak and Kodis 1983). A more modern version of RADFLOz is embodied within 
the code HYCHEM, which includes the nonequilibrium chemistry. RADFLOz and 
HYCHEM have been modified and improved extensively during the past year. 
Further code improvements and extensions in the 1980s included the development 
of HYCHEM. HYCHEM is an extension of RADFLO that includes computations of gamma 
ray and neutron-produced ionization rates with associated nonequilibrium chemistry 
processes and uses the computed chemical species concentrations in computing the radiative 
absorption coefficients, thus fully coupling the chemistry and radiation transport 
computations. As a test of HYCHEM, the 1956 Redwing-Navajo chord experiment was 
modeled, including time-resolved spectroscopic measurements of the gamma ray and neutron- 
produced smog outside of a fireball (Zinn et al. 1980, 1982). Several related studies validated 
the chemistry. These included development of a Teller light model in the TELLER and 
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SLEMP computer codes, which were used to model the time-resolved Teller light- 
spectroscopy experiments in the Redwing-Dakota event (Zinn and Sutherland 1987). 
Other code modifications in the late 1980s included development of a version with 
a noniterative Godunov-Riemann equation-solver hydrodynamics subroutine in place of 
the explicit Richtmyer-von Neumann hydrodynamics algorithm of previous versions (von 
Neumann and Richtmyer 1950). Also, a multimaterial RADFLO version (Symbalisty 1994) 
treats the bomb materials and possible surrounding solid or liquid layers as distinct from 
the surrounding air (i.e. with distinct equations of state and opacities). In addition, a research 
study was undertaken to incorporate effects of Taylor instabilities in the bomb debris 
expansion (Jones and Kodis 1986; Jones et al. 1979, 1981) but was discontinued. 
Some further history relating to equations of state and absorption coefficients is 
included in section 3.1. 


3.0 THE MODEL 
We list here the radiation-hydrodynamic equations that the RADFLO code solves 
(Zinn 1973) in spherical symmetry. The equations of conservation of mass, momentum and 
energy are 


∂ρ 
∂τ = −∇ ⋅ (ρv) 


∂(ρvα ) 
∂t 
= −∇⋅(ρvαv) − ∇αP 


∂e 
∂t = −∇⋅(ev) − P∇⋅ v − 
0 


∞∫ ∇⋅ Fvdv 


where the monochromatic flux function is 


Fv = 2π 
0 


π∫ Iv(θ)sinθ cosθdθ 
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and the specific monochromatic intensity at frequency v is 


Iv(P2) = Iv(P1)exp − 
′ 
µ (ξ)dξ 
P1 


P2∫ 


+ 
′ 
µ (ξ)Bv 
P1 


P2∫ 
(ξ)exp − 
′ 
µ ( ′ 
ξ )d ′ 
ξ 
ξ 


P2∫ 
dξ . 


P1and P2 are two points on the surface of a volume element. The angle θ is the 


angle between the radius vector and the direction of a specified ray. The differentials dξ 


and d ′ 
ξ are elements of distance along the ray connecting points P1and P2. Bv is the 
Planck function and ′ 
µ is the linear absorption coefficient corrected for stimulated emission. 
The pressure, P , and temperature, T, are assumed to be known functions of the 
mass density, ρ , and the material internal energy density,e: 


P = P(ρ,e) 
, 
T = T(ρ,e). 


The radiative absorption coefficients are also assumed to be known: 


′ 
µ = 
′ 
µ (ρ,e,v). 


Note that the force of gravity is not included in these equations, so that buoyant deformation 
of the fireball is not considered (nor are other two-dimensional effects). This implies a limit 
to the interval of time for reliable results. However, for low-altitude explosions this interval 
generally encompasses the time during which the main fireball radiation is emitted, unless 
there are significant ground interactions. 
The basic radiation transport algorithm used in RADFLO is unchanged from that 
described in (Zinn 1973), except for two details which will be described below. The 
Richtmyer-von Neumann hydrodynamics algorithm is also unchanged. However, we have 
added the Godunov hydrodynamics subroutine which we now use as the standard. The older 
Richtmyer-von Neumann hydro-subroutine is still retained in the code for optional use. The 
set of air equation of state and opacity tables is unchanged except for a few minor changes 
in the opacities. The bomb initialization algorithm is also basically the same, although we 
now allow for an arbitrary number of mesh cells within the bomb volume. A large amount 
of new code was added in 1981 for patching in the SSW results, and these sections have 
been modified further during 1994–1995. 
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3.1 Air Equation of State and Absorption Coefficients 
The RADFLO, HYCHEM, and SSW codes all contain built-in tables of Planck 
functions and air equations of state and opacities. (The opacity 
′ 
κ [cm2/g] is equal to the 
linear absorption coefficient 
′ 
µ (cm-1) divided by the density ρ —both at specified frequency 


v and corrected for stimulated emission). 
The equation of state tables are used to find the cell values of the pressure P and 
temperature T from the corresponding computed cell values of ρ and E. For ease and 
accuracy of interpolation the equation of state data are stored in the form of P / ρE and 


T / E vs E andρ , using evenly-spaced logarithmic arrays of E (90 values from 2 ×109 to 


1×1016 erg/g) and of ρ (eight values from 1.29 ×10−9 to 1.29 ×10−2 g/cm3). Since the 
logarithmic spacings in both arrays are fixed, the interpolations can be very rapid, since for 
given values of a local cell energy and density the interpolation parameters can be calculated 
directly, without table searching. The quantities P / ρE and T / E are relatively slowly 
varying functions of ρ and E, so that interpolation of those quantities instead of just P 


and T can be carried out with relatively little error. When P / ρE and T / E have been 
interpolated, the corresponding P and T can be immediately calculated. 
The opacities are stored in a three-dimensional array as functions of ρ , T, and hv. 
The same evenly-spaced logarithmic array of ρ is used, together with an evenly-spaced log 
array of T values and an evenly-spaced array of hv values. Interpolation of opacities, ′ 
κ , 
instead of 
′ 
µ minimizes interpolation errors, since 
′ 
κ varies more gradually with density 
than does ′ 
µ , and as in the case of the equation of state data above, the interpolation indices 
can be calculated directly without table searches. With respect to the quantum energy variable 
hv, no interpolations are needed since the radiation transport computation is carried out for 
each of the same 42 groups as those of the opacity tabulation. Once the interpolated values 
of ′ 
κ are obtained, the absorption coefficients can be immediately calculated from ′ 
µ = ρ ′ 
κ . 
For each of the logarithmically-spaced values of hv in the tabulations there 
corresponds a quantum energy interval δhv. For each T and hv we also tabulate a 


precomputed value of the Planck function integral 


hv1i 


hv2i 
∫ 
Bvdv, where hv1i and hv2i are the 


values of hv at the two ends of the interval i. In carrying out the transport computation, 
when Planck function integrals are required we interpolate them from the tables in the same 
manner as the opacities. This is faster computationally than computing the Planck integrals 
on the fly. 
The same interpolation scheme and the same tabular arrays of ρ , E, and hv have 
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been used in all versions of RADFLO and HYCHEM since 1965. Also the same equation of 
state tables are used. However, between the RADFLOz and RADFLOs versions there are 
some small differences in the opacity tables. 
With respect to development of the opacity data, gaps exist in the history due to 
poor documentation. Intensive research efforts to produce high-temperature air equation of 
state data (and computations) began during the Manhattan Project (Fuchs et al. 1942), and 
these efforts increased through the mid-1960s. In the 1950s and 1960s, most of this work 
was funded by the U.S. Air Force through the Air Force Special Weapons Command 
(AFSWC) and later the Air Force Weapons Laboratory (AFWL) and was carried out primarily 
at the RAND Corporation (Gilmore and Latter 1956; Gilmore 1955 and 1958; Latter 1954), 
and the National Bureau of Standards (Hilsenrath, Green, and Beckett 1956; Hilsenrath, 
Klein, and Woolley 1959; Woolley 1962). The equation of state was fairly well defined by 
the time of the initial development of RADFLO in 1963. 
Research to produce reliable air-opacity data extended from the 1940s through the 
mid 1970s. Considerable work to assemble opacity tables for use in a fireball radiation 
transport code was done by A. Skumanich (1960) with help from H.A. Bethe. A large part 
of the opacity data came from Gilmore and Latter (1956). The first RADFLO code borrowed 
heavily from Skumanich’s work. The first RADFLO opacity tables included only continuum 
absorption processes including free-free and bound-free transitions from atoms and ions of 
nitrogen, oxygen, and argon. No molecular bands or atomic lines were included. Below 
the atomic oxygen L edge, the bound-free absorption coefficients for hot air were estimated 
using the Raizer approximation and with advice from H.A. Bethe. 
From the 1950s through the 1970s, research on opacities was supported by AFSWC, 
AFWL, Los Alamos Scientific Laboratory, the Air Force Office of Scientific Research, Naval 
Research Laboratory, Defense Atomic Support Agency, and DNA. Experimental and 
theoretical work was carried out at RAND Corporation (Gilmore and Latter 1956; Gilmore 
1958), Lockheed (Churchill, Hagstrom, and Landshoff 1964; Churchill et al. 1966; 
Buttrey and McChesney 1965; Mueller 1963; Johnston and Platas 1969; Armstrong et al. 
1958; Churchill and Hagstrom 1964; Churchill et al. 1962; Johnston et al. 1971), AVCO 
Corporation (Allen et al. 1962; Keck et al. 1963; Taylor 1960; Taylor and Caledonia 1968; 
Wray and Taylor 1969; Keck et al. 1959; Kivel and Bailey 1957; Camm et al. 1967; Wentink 
et al. 1967;), Los Alamos Scientific Laboratory (Mayer 1947; Cox and Stewart 1964), AFWL 
(Harris and Sullo 1970; Generosa and Harris 1970; and Harris et al. 1969), General Electric 


11 


Corporation (Breene 1966), California Institute of Technology (Patch et al. 1962), and 
elsewhere. Funding for opacity work tapered off sharply in the 1970s. 
For the original RADFLO code, which included only continuum contributors to the 
opacity, 42 frequency groups were considered to be sufficient to resolve the most important 
structure in the spectra. However, as time went on, it became clear that accuracy would be 
improved greatly if molecular band systems and atomic lines were included. The fine spectral 
structures of these band systems and lines could not be resolved with 42 frequency groups 
so it was necessary to decide how to average the spectral details over each of the frequency 
groups. For spectral regions, temperatures, densities, and mesh cell thicknesses for which a 
cell was optically thick, a reasonable procedure was to use a “Rosseland” average of the 
absorption coefficients over each of the frequency groups. The Rosseland average tends to 
weight the valleys between the spectral lines. On the other hand, for conditions in which a 
cell is optically thin and where emission is dominant over absorption, generally it was better 
to use a “Planck” mean, which approximates a direct average of the absorption coefficient 
over the given frequency group. 
In the late 1960s opacity groups were formed in the Los Alamos Theoretical 
Division and at AFWL. They maintained a close collaboration. The AFWL group produced 
a computer code known as ABSCO, which assembles atomic and molecular line atlases for 
given conditions of temperature and density (thermodynamic equilibrium assumed). The 
Los Alamos group produced a code that forms specified types of spectral averages of the 
opacities over specified frequency groups. This provided the opacity tables that we currently 
use in RADFLOz (and HYCHEM). These opacities are Rosseland averages over each of 
the specified 42 frequency groups that we use. The atomic, ionic, and molecular processes 
included in this opacity set are 


• 
Free-free absorption in the fields of ions and neutrals, given the prevailing equilibrium 
populations of each. 
• 
Bound-free continua for N atoms and ions, corresponding to the prevailing ionization 
levels and quantum state populations. 
• 
Bound-free continua for O atoms and ions, corresponding to the prevailing ionization 
levels and quantum state populations. 
• 
Bound-free continua for Argon atoms and ions, corresponding to the prevailing 
ionization levels and quantum state populations. 
• 
Bound-bound transitions of N and O. 
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• 
Photodetachment from O– and O2 
–. 
• 
Photoionization of N2 and O2. 
• 
O2 Schumann-Runge continuum. 
• 
O2 photodissociation continuum number 2. 
• 
O2 Schumann-Runge band system. 


• 
N2 
+ band systems (9). 
• 
N2 band systems (4). 
• 
NO band systems (11). 
• 
NO+ band system (1). 
• 
CO band system (1). 
• 
NO2 bands. 
• 
CO2 UV band. 
• 
CN band systems (3). 


Except for addition of the CO absorption the RADFLO opacities have not been 
modified since 1974. 
Minor differences in the opacities used in RADFLOs are described in Sowle et al. 
(1972) and Sappenfield (1973). In particular, the RADFLOs opacities are 30% to 50% smaller 
than the RADFLOz opacities in the temperature-wavelength regime that determines the 


optical power minimum time (i.e. kT = 0.2 to 0.6 eV and hv = 1.5 to 4 eV). In this range the 
absorption is mainly due to NO2. 
In Figures 1–4 we plot the Rosseland and Planck mean absorption coefficients for 
eight different air temperatures as a function of photon energy. In RADFLOs, the Rosseland 
mean is used if the optical depth, computed with the Rosseland values, is greater than one. 
The Planck mean is used if the optical depth, computed with the Planck values, is less than 
one, and a geometric average is used otherwise. The Planck mean absorption coefficient is 
always greater than or equal to the Rosseland mean. 
The RADFLOs equation of state (EOS) for air is the MRC equation of state for air 
developed by Sappenfield and is very similar to the RADFLOz EOS. It relates the density 
(ρ ), the internal energy (e), the temperature (T), and pressure (P). It is based on equation of 
state tables of Hilsenrath et al. (1956, 1959) and Gilmore (1955). The temperature range 


of validity of the tables is 1000 K to 5 ×106K. Extrapolation beyond this range is 


straightforward. In terms of density the tables are valid for 1.29 ×10−9 ≤ ρ ≤ 1.29 ×10−2 gm/ 
cc. 
The EOS will return (T,γ −1)when given the state variables (ρ ,e ) (Figures 5 and 
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6). The pressure is then p = (γ −1)ρe . The EOS will also return (γ −1,e) when given 
(T,ρ). 
Godunov’s method in the hydrodynamic solver requires two more parameters derived 
from ideal gas relations and the equation of state, the sound speed, cs, and a shock parameter, 


ra. These parameters are defined as follows: 


cs 
2 = γP / ρ, and ra = 1+ (γ −1) / 2. 


3.2 Godunov Hydro 
The staggered-mesh, Richtmyer-von Neumann hydrodynamics algorithm has been 
replaced with a variation of the Godunov method and uses an approximate Riemann solver 
(Dukowicz 1985). In fact, we are using a one-dimensional, spherically symmetric analog of 
the two-dimensional CAVEAT code developed at Los Alamos (Addessio et al. 1992). We 
quote from Addessio et al. (p. 7): 


The benefits of a Godunov approach are of considerable importance. Staggered- 
mesh techniques result in a close coupling between pressure and velocity fields but 
usually require artificial dissipation for numerical stability and/or the achievement 
of proper entropy changes. Such artifical dissipation may be explicit (artifical 
viscosity) or hidden (donor-cell advection), but such methods typically fail 
to overcome mesh-drifting tendencies and are often excessively diffusive. With the 
Godunov approach, cell-edge pressures and velocities required for acceleration and 
work terms are obtained by solving the Riemann problem each cycle from the cell- 
centered state variables on the two sides of each cell edge. The result is again close 
coupling between pressure and velocity fields, usually with the optimal local level of 
dissipation. Because the minimal amount of dissipation is added by the Riemann 
solution, the Godunov method results in lower overall diffusion, and this does not 
mask the effect of realistic sources of diffusion from physical processes. 
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Figure 1. The Rosseland mean absorption coefficient as a function of photon 
energy for four different temperatures from the RADFLOS tables. The air tempera- 
tures, in eV, are 0.862, 1.034, 1.206, and 1.379. 
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Figure 2. The Planck mean absorption coefficient as a function of photon energy for 
four different temperatures from the RADFLOS tables. The air temperatures, in eV, 
are 0.862, 1.034, 1.206, and 1.379, which are the same values as the previous figure. 
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Figure 3. The Rosseland mean absorption coefficient as a function of photon energy 
for four high temperatures from the RADFLOS tables. The air temperatures, in eV, 
are 412, 498, 600, and 725. 
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Figure 4. The Planck mean absorption coefficient as a function of photon energy for 
four high temperatures from the RADFLOS tables. The air temperatures, in eV, are 
412, 498, 600, and 725, which are the same values as the previous figure. 
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1.20E–03g/cc, 412.0, 498.0, 600.0, 725.0 eV 
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Figure 5. Plots of the internal energy vs temperature for eight atmospheric densities. 
The densities, left to right, are 10-2, 10-3, 10-4, 10-5, 10-6, 10-7, 10-8, and 10-9 gm/cc. 
Ambient sea level density is ≈ 1.23 ×10−3gm / cc. 
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Figure 6. Plots of (γ -1) vs temperature for eight atmospheric densities. The densi- 
ties, top to bottom, are 10-2, 10-3, 10-4, 10-5, 10-6, 10-7, 10-8, 10-9 gm/cc. 
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We next describe our one-dimensional analog of CAVEAT. CAVEAT is a time- 
explicit, finite volume formulation of the Navier-Stokes equations. We ignore viscosity, heat 
conduction, and body force (gravity) terms. The differential equations, ignoring the radiative 
flux terms, are integrated over a finite volume element defining a discretized mesh. In one- 
dimensional, spherical symmetry the volume elements are spherical shells. Let 


1. UC = the cell centered velocity 
2. r = the cell face radius 
3. E = the cell centered total energy 
4. e = the cell centered internal energy 
5. P = the cell centered pressure 


6. V = the cell volume = 4π(ri+1 
3 − ri 
3) / 3 


7. A = the cell face area = 4πr2 


8. CM = the cell mass = constant for each cell in a Lagrangian simulation 
9. DR = the cell width = ri+1 − ri 


The approximate Riemann solver—given the cell centered state variables, cell geometry, 
and an equation of state—computes 


10. UF = the cell face velocity 
11. PF= the cell face pressure 


which are used to update, to time tn+1, the state variables of cell i by 


ri 
n+1 = ri + ∆tUFi 


UCi 
n+1 = UCi − ∆t(PFi+1 − PFi)Vi / (CMiDRi) 


Ei 
n+1 = Ei − ∆t((PF ⋅ A⋅UF)i+1 − (PF ⋅ A⋅UF)i) 


where everything on the right-hand side of the three above equations is evaluated at . ∆t is 


the computational time step equal to tn+1 − tn. The new cell geometry can now be calculated 
and the mass density and the internal energy follows: 


ρ = CM / V , e = E / CM − UC2 / 2. 
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3.3 The SSW Patch 
The much greater complexity of the present RADFLO codes compared to the 1973 
version is mainly associated with the “SSW patch”—the extra logical switches and 
computations necessary for deciding where the fireball front is, what its velocity is, how to 
correct for artificial viscous pressure or other numerical diffusive effects, and how to correct 
the computed radiation fluxes out of the fireball surface to correspond to the SSW values. It 
is beyond the scope of this report to describe the detailed reasons for the numerical oscillations 
in brightness and radiated power which the SSW patch is designed to avoid. They are 
discussed in detail in Zinn and Sutherland (1982). In brief outline they arise because of 
numerical smearing of the fireball front, which occurs in two different ways: 


1. Relative to the dimensions of the numerical mesh cells (meters) the early fireball “front” 
is essentially a step function in temperature and density. The hot air behind the front is 
opaque to radiation at visible and IR wavelengths. In the computation, as the model 
front progresses through the mesh it is necessarily averaged at each instant across one 
or more cells, so that the cell that contains the front has an artificial temperature that 
increases with time from the initial ambient temperature up to the true fireball front 
temperature. 


2. The hydrodynamics algorithm intentionally smears the shock front over several cells 
since it can’t cope with a true discontinuity. In both the Richtmyer-von Neumann and 
Godunov hydro algorithms, the artificial smearing processes lead to a poorly determined 
fireball surface temperature and to oscillations in computed brightness. 


To cope with the problem we had to devise ways to determine (1) where the actual fireball 
front is and (2) what its true temperature is, or alternatively, what its velocity is. 
Both the artificial viscous pressure construct and the Godunov method have 
the important property that they produce the correct jump in density and entropy (or 
temperature) across the model shock front. They also produce a smooth transition in energy 
across the front, instead of a discontinuity. In doing so, they produce an artificially elevated 
energy and temperature just ahead of the shock— where is positive. The shock “front” lies 
within a cell where ∂ρ / ∂t > 0 but where ∂ρ / ∂t ≤ 0 in the cell just behind. In the cells just 
ahead of the shock (where ∂ρ / ∂t > 0) there is artificial compressional heating, and the 
temperature is artificially elevated, while behind the shock (where ∂ρ / ∂t ≤ 0) the 
temperature is correct. Thus, after each hydrodynamics step of the computation and 
before each radiation transport step, we must correct the cell internal energies ( E) for the 
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artificial heating. To do this in the code, we maintain two separate cellwise energy arrays, the 
first being E , the energy used by the hydrodynamics subroutine, and the second being 


eshk , the energy used by the radiation transport routine. E is updated each time step, both 
by the hydrodynamics and radiation transport routines. eshk is identical to E in the cells 
where ∂ρ / ∂t ≤ 0, but in cells where ∂ρ / ∂t > 0 it is updated only by the radiation transport. 


eshk is a better approximation to the true internal energy than is E (since it does not contain 
the artificial viscous heating). Still, as the necessary consequence of using a finite mesh, the 
transition in eshk at the fireball front is spread artificially over two or more cells. The 
temperatures used by the radiation transport routine are computed from eshk and the density 
rho—these temperatures are used in interpolating the absorption coefficients and computing 
Planck functions. For purposes of computing the fireball growth and profiles of temperature, 
density, velocity, etc., no other corrections are made for the numerical smearing of the front. 
However, for computing the fireball radiated power and visible and thermal brightnesses it 
is necessary to go a step further. 
Because of the spreading of the fireball front temperature discontinuity over two or 
more cells, the brightness and power radiated by the front will still tend to oscillate as the 
front progresses through the mesh. The oscillation goes through one cycle for each mesh 
cell traversed. In HYCHEM (or RADFLOz) to deal with this problem and to compute the 
“correct” brightnesses for wavelengths longer than 180 nm, we proceed as follows: as the 
computation runs we attempt to determine at each numerical time step (1) where the true 
fireball front is and (2) what its velocity is. From this velocity and a stored tabulation of 
effective radiating temperature vs front velocity (as derived from the separate SSW 
computations) we interpolate to obtain the effective radiating temperature of the front. This 
radiating temperature includes the absorption by the fine-scale radiative precursor, which 
is resolved in the SSW computations but is much too thin to be resolved in the coarse zones 
used by RADFLO. This effective radiating temperature is used to compute the 
effective blackbody flux out of the fireball surface for each of the RADFLO quantum energy 
groups (but limited to wavelengths longer than 180 nm). The radiation transport routine is 
run in the usual way, except that we substitute the SSW effective blackbody flux at the outer 
vertex of the cell that contains the front. This flux may then be modified by the normal 
transport algorithm to account for absorption by preheated air outside the fireball front (i.e. 
the x-ray veil, and/or by “smog” in the case of HYCHEM computations). As a further fussy 
detail, we also modify the computed absorption that arises from the cell containing the front 
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because that cell has a distorted temperature arising from the abovementioned front-spreading. 
Our previously published SSW computations were only for shocks and thermal waves 
in sea level-density air. However, in order to use the SSW patch algorithm in HYCHEM 
computations for altitudes above sea level, we have recently run additional SSW 
computations for a range of air densities. The results of these computations are outlined in 
Appendix A. Qualitatively they show that as the air density decreases the shock velocity for 
maximum shock brightness decreases, and the magnitude of the maximum brightness 
also decreases. On the other hand, in the high-velocity thermal-wave regime the wavefront 
brightness for a given velocity increases with decreasing air density. These results are 
incorporated in tables within HYCHEM, and the SSW patch algorithm has been extended 
to include the background air density as a variable in addition to the wavefront velocity. 


3.3.1 The RADFLOs Patch 
The RADFLOs code also contains a fireball brightness patch algorithm which is 
designed to give an approximate fit to the SSW results from Zinn and Sutherland (1981), 
relating to the fireball front radiative precursor attenuation for sea level air-density cases, 
but at the same time provide a guess as to the amount of precursor attenuation that would be 
expected in high-altitude cases. That algorithm of course predates our recent extended set 
of SSW computations which include other than sea level air densities. 


3.4 Time Steps 
The basic numerical time step dt is dictated by requirements of numerical stability, 
as described by the familiar Courant conditions for the hydrodynamics and analogous 
conditions for the radiation transport. Generally it takes many time steps for the fireball 
front to progress from one cell to the next. We also refer to a “macro” time step—namely 
the time required for the front to cross a cell. dt refers to the “micro” time step. At each 
micro time step we determine the index lzero of the cell that contains the fireball front. The 
algorithm for defining lzero is fairly complex, but basically it is the larger of two indices, 


irmx and itmx , where irmx labels the outermost cell with a density that exceeds a certain 
threshold indicating an “air shock,” and itmx labels the cell that contains the largest 
temperature gradient. As the model run progresses, the times and radial positions at which 


lzero advances from each cell to the next are stored in a pair of arrays, and the fireball 
front velocity is obtained as the finite difference derivative—one value for each macro time 
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step. These velocities are used to interpolate in the SSW-generated tables of effective radiating 
temperature vs velocity to obtain the effective fireball temperature for the brightness 
and radiated power computation. The radiated flux patch is applied at the cell vertex just 
outside of lzero. 


3.5 Chemistry 
In HYCHEM, if the input parameter ichem is set to 1, the run will include 
chemistry—a set of chemical-reaction kinetics equations is integrated for a predetermined 
set of chemical species (currently 45) for a predetermined set of time steps. The results 
include profiles of the concentrations of these species at the predetermined output times. The 
chemistry outside the fireball is driven largely by neutron- and gamma ray-induced ionization 
of the air and ensuing chemical reactions. Therefore the code necessarily includes algorithms 
for computing the prompt neutron and gamma ray emission and transport and the bomb- 
debris radioactive-decay gammas. An important effect of the chemistry is to produce 
optically-absorbing chemical species outside the fireball, which tend to absorb radiation 
emitted by the fireball and thereby reduce the radiated power. These species include NO2, O3, 
and HNO2, and in fireball vernacular they are called “smog.” 
The chemistry computations are done at prescribed intervals of time which are much 
larger than the micro time step used in the main rad-hydro computation. After the chemical 
concentrations have been computed, we can compute the changes in the radiative absorption 
coefficients resulting from the deviations from chemical equilibrium. The equilibrium 
chemical composition is determined by a separate subroutine that computes the composition 
that gives the minimum value of the Gibbs free energy. Then at each mesh cell at each micro 
time step, the difference between each concentration and the equilibrium value is used 
together with a set of tables of absorption cross sections to correct the normal 
equilibrium radiative-absorption coefficients. In this way the radiation transport is modified 
by the chemistry, and the effects of nuclear smog outside the fireball are computed. The 
smog absorption can have a very large effect on the radiated power from a large-yield fireball 
at early times. 
The chemical reaction rate equations include photochemical reactions associated 
with radiation emitted by the fireball. The reaction rates are recomputed at each chemistry 
time step to take account of temperature changes. The chemistry and gamma ray+neutron 
transport algorithms are described in some detail in Zinn et al. (1982) along with the chemical 
species set and the reactions and rate coefficients. The algorithm for the neutron transport in 
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the present HYCHEM code is quite crude and probably constitutes the most serious source 
of error in the chemistry and smog computations. At one time we had hoped to add a detailed 
neutron transport subroutine based on the standard code MCNP. 


3.6 The RADFLOs Smog Model 
The RADFLOs code also includes an algorithm for smog absorption, although no 
chemistry computations are done. This algorithm (Sappenfield 1976) is based on the analysis 
of data by C. K. Mitchell from the Navajo and Flathead chord experiments, in which it was 
assumed that the smog optical-absorption coefficient was directly proportional to the time- 
integrated gamma ray and neutron energy deposition, with no other time dependence 
(Mitchell 1979). The gamma ray- and neutron-energy deposition is computed with an 
algorithm from the DNA ROSCOE code (Knapp et al. 1974). 


3.7 Debris Opacities 
Since the tabulated opacity data are for air, the radiation transport in the early bomb 
vapor and later bomb-debris-air mixture are treated only approximately. During the early 
period of x-ray emission by the bomb, we multiply the absorption coefficients within the 
debris by an arbitrary factor of ten to assure that the bomb emits as a blackbody. For later 
times we add extra terms to the opacity in the debris region to account for the expected 
higher free electron concentrations (due to the lower ionization potentials of debris materials) 
and the expected large number of allowed low-energy bound-bound transitions. For the 


bound-bound transitions we assume an average absorption cross section of 1.3 ×10−18cm2 


over the range hv = 1 to 4 eV, which assumes that each debris atom has one allowed transition 


in this range with an oscillator strength of 107s−1. The extra terms for the enhanced free 
electron concentration include (1) enhanced free-free absorption in the fields of neutral atoms 
and (2) enhanced bound-free absorption due to negative O- ions. We also estimate the size 
of the debris region, since the debris expansion is assuredly Taylor-unstable, using an 
algorithm that says that the debris tends to mix with air up to 20 times its own mass. 


3.8 Accessory Subroutines 
In addition to the main computations of fireball growth and chemistry, the code 
includes accessory routines for things such as (1) computing profiles of fireball brightness 
at visible wavelengths, (2) storing tables of fireball radii, velocities, and radiated power vs 
time for later plotting, and (3) input and output. The code also includes subroutines for 
adjusting the mesh as the problem runs—i.e., (1) combining cells behind the front to improve 
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running speed, and (2) adding cells at the outside of the mesh when the shock wave gets 
close to the outer boundary. These routines are necessary to enhance running speed 
and spatial/temporal resolution. 


4.0 INITIAL CONDITIONS 
The yield of a given device is nominally partitioned in terms of an x-ray yield, a 
kinetic energy (of the debris) yield, and an internal energy (of the debris) yield.1 In addition 
a gamma ray yield and a neutron yield may be specified. 
In RADFLOz, the algorithm for initialization of the explosion is based on a very 
simple model for a generic nuclear bomb and has been described previously in Zinn (1973). 
We specify the device yield (Y), the mass (M), the average density (ρs), and a time ( τ x) 
representing the time duration of the nuclear reaction. We further assume that the device is 
spherical, and that the energy (Y) is deposited uniformly within the device mass so that the 
density and temperature remain uniform. The x-rays from the explosion are emitted at a 
constant rate over the time ( τ x ) with a Planckian distribution characterized by a 
fixed temperature Tx (which is to be computed). In addition we assume that the 
energy remaining after x-ray emission is apportioned equally between internal and kinetic. 
To compute Tx we also need to know the effective specific heat for raising the bomb mass 
to temperature Tx , and the radiating surface area, which is obtained from the given mass 
and given density. Given these assumptions it is then straightforward to calculate the 
radiating temperature Tx and the partition of the energy Y between (1) radiated, (2) internal, 
and (3) kinetic. The equations are given in Zinn (1973). 


4.1 X-ray Deposition 
A special subroutine is used for computing the distribution of energy deposition in 
the air due to the bomb x-rays. In the case of a nuclear explosion at sea level, the electron 
collisions are fast enough to assure that the air stays in thermodynamic equilibrium as the x- 
rays are deposited. The air-absorption coefficients change as the air heats up, but they can 
be computed in the usual way from the equilibrium air-opacity tabulation. The evolution of 


1 The word debris is synonymous with device or bomb in this report. In RADFLOs, the first jdeb number of 
zones are debris zones. The user must supply an internal energy and kinetic energy to each of them. The 
default equation of state for the debris is the same as the air equation of state. 
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the deposition process must be done stepwise to follow the dynamically changing 
absorption coefficients. At each step the x-ray energy that is transported and deposited in 
the air is, at the same time, removed from the bomb. The radiation transport is computed in 
the usual way by the algorithms already described. The stepwise computation is continued 
until the previously determined x-ray yield has been transferred from the bomb to the air. In 
the case of a high-altitude explosion (about 40 km or higher), electron collisions are not 
sufficiently fast to maintain thermodynamic equilibrium on the time scale of the bomb x- 
ray emission. For such cases it is a better approximation to assume that the air-absorption 
coefficients are those of cold air—except where the energy deposition density is so large as 
to exceed the energy of four x-ray photons per air atom, in which case the absorption 
coefficient drops to a very small value (the bremsstrahlen limit). The four-photon limit is 
based on the expectation that each x-ray photon absorbed (in the K shell) results in ejection 
of two electrons—a photolectron and an Auger electron, and the K shell is immediately 
refilled from outer shells until there are no more outer shell electons. In the case of oxygen, 
which has eight electrons per atom, this permits four photons to be absorbed. 


4.2 HYCHEM Input Parameters 
In the input to the HYCHEM code the bomb, or device, is described by 12 parameters 
as listed in Table I. The most important of these are (1) the yield (yield), (2) the bomb mass 
(bmbms), (3) the average bomb density (rho1), (4) the x-ray radiating time (taux), and (5) 
the fraction of the non-x-ray yield (or non-radiated) yield that ends up as kinetic energy 
(cekin/2 ). In addition there are three prompt gamma ray parameters, namely (6) the fraction 
of the device yield emitted as prompt gammas (fgam), (7) the approximate gamma ray straggle 


time (dtgam), and (8) the average gamma ray mass absorption coefficient (akappa). Also, 
there are three prompt neutron parameters, namely (9) the neutron yield fraction (fneut), (10) 
the average neutron velocity (vneut), and (11) the average neutron mean free path 
(alam). Finally, there is (12) the fission yield fraction (ffiss), which is used for computing the 
delayed fission-product decay, beta- and gamma-emission rates. If chemistry is not included 
in the computation, only the first five of these parameters are relevant. In the input subroutine 
all the bomb parameters except yield are provided with default values, in case they are not 
otherwise specified. The default values are also listed in Table I. 
The computational mesh is set up with an automatic algorithm which is controlled 
by four input parameters, namely: (1) the number of cells containing bomb material (ndc), 
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Table I. HYCHEM Input Variables. 


Variable 
Meaning 
Default Value 


akappa 
The average gamma ray mass absorption 
coefficient in cm2/gm. 
0.0113 


alam 
The average neutron mean free path in cm. 


bmbms 
The bomb mass in gm. 


cekin 
The fraction of the non x-ray (or non 
radiating) yield that ends up as kinetic energy 
is cekin/2. 


1.0 


dtgam 
The approximate gamma ray straggle time in 
s. 


dtk 
Time step limiter parameter. 
0.2 


eint2 
Initial ambient atmospheric internal energy in 
erg/gm. 


endtime 
Final problem time in seconds. 
computed if 0 


ffis 
Fission yield fraction. 
1.0 


fgam 
The fraction of the device yield emitted as 
prompt gammas. 
0.001 


fneut 
The neutron yield fraction. 
0.01 


i64bit 
32 vs 64 bit machine switch. 
0 


ichem 
Chemistry switch. 
0 


ifork 
Restart switch. 
0 


ihialt 
High-altitude switch. Zero for low low-altitude 
explosions and 1 for above 40 km. 
0 


ihist 
Chemistry diagnostics switch. 
0 


imax 
Maximum number of cells. 


iriem 
Riemann hydro switch. 
0 


isswfit 
Switch for SSW patch. 
1 


ndc 
Number of debris cells for the problem. 


nq 
Viscosity coefficient. 
2 


nzone2 
Number of zones in the fine zoned region after 
the debris cells. 


qb 
Viscosity coefficient. 
1.0 


rho1 
The average bomb density in gm/cc. 
2.0 


rho2 
Initial density of air cells in gm/cc. 
0.0012 


rmax 
Maximum initial radial extent of the mesh. 
computed if 0 


taux 
The x-ray radiating time in s. 


tchm 
Time at which chemistry is turned on. 


tchmfrc 
Fraction of the elapsed time before doing 
chemistry again. 


vneut 
The average neutron velocity in cm/s. 


xh2o 
Mole fraction of water vapor, used if 
chemistry is turned on. 


yield 
The yield in kilotons. 
none 


4 ×104 


105 yield1/3 


10−7 


10−8 


5.2 ×109 
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(2) the number of cells between the bomb surface and the radius at which the average fireball 
temperature is about 30 eV (nzone2), (3) the total number of cells out to the outermost mesh 
boundary (imax), and (4) the radius of the outermost boundary (rmax). If rmax is not otherwise 


specified it is set with a default expression r max = 780 × (yield / rho2)1/3cm2. rho2 is the 
background air density (g/cc), which must also be specified. The 30 eV radius is computed 


as r30eV = 69.3 × (yield / rho2)1/3. 
The zoning algorithm computes the initial cell vertex radii so that (1) the ndc bomb 
debris cells contain equal masses, whose sum is bmbms; (2) for the next nzone2 cells, each 
cell mass is a constant multiple dmfctr2 times the mass in the previous cell; (3) likewise, for 
the remaining nzone3 cells, each cell mass is a constant factor dmfctr3 times the mass of the 
previous cell. 


4.3 RADFLOs Input Parameters 
RADFLOs assumes the user of the code is fairly knowledgeable about the device 
being modelled and thus has a large number of input parameters for fine tuning the simulation. 
We list them in Tables II and III and will describe them in detail in the User's Manual report 
(Horak 1983). 


5.0 SOME NUMERICAL RESULTS 
In order to test the current code and to benchmark a reference set of computed results, 
we have run a number of cases. Most of them were run with HYCHEM (or RADFLOz), but 
a few were run with RADFLOs. The cases include yields ranging from 0.01 kilotons to 1 
Megaton, and altitudes ranging from sea level up to 50 km at a few different device masses. 
In most cases we have not included the nonequilibrium chemistry (standard RADFLO mode), 
but in a few instances we ran identical cases both with and without the chemistry. 
Since this is an unclassified report which describes the computer code and general 
aspects of fireball physics, we will avoid any quantitative references to nuclear test data. 
The computations are all based on our hypothetical model of a generic bomb, as described 
above. In most cases we have used the default values of bomb mass, etc. However, in a few 
cases we investigated effects of varied bomb masses. In RADFLOs cases we attempted to 
match the input parameters as closely as possible with an equivalent run done with RADFLOz, 
so that the results would be directly comparable. 
In this section we describe results relating to optical observables—especially optical 
power vs time (bhangmeter curves), and fireball radius vs time. 
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Table II. RADFLO Input Variables. Most of the variables are defined as they were 
in Sappenfield (1984). 


Variable 
Meaning 
Default Value 


cno2 
Factor for determining 
concentrations in 
smog routines. 
Depends on 
fgfact/fnfact 


dfrac 
The fraction of the shock radius at which the 
outer boundary of the debris is assumed to be, 
prior to the time at which the debris boundary is 
frozen in Langrangian space. 


1.0 


dr 
Gridding array. 


dth 
Initial hydro time step. 


ekin0 
Initial kinetic energy in kilotons. 
none 


eofx(jdeb) 
Initial bomb internal energies. 


est 
Ergs/gm required to strip an air atom at 
equilibrium. 


fact 
Total gamm ray and neutron yield (ergs). 
0.0 


fcng1 
The maximum fractional internal energy change 
per time step due to radiation flow that is 
allowed initially. 


fcng2 
The value to which fcng1 is changed at tfcng1. 


ffrc 
Fission fraction. 
0.5 


fgfact 
Fractional yield in gamma rays. 
0.0 


fh 
Fraction of total energy in the spectrum with 
temperature th. 


fnfact 
Fractional yield in neutrons. 
0.0 


i 
Flag for new start or restart. 
none 


idb 
Flag indicating treatment of debris contribution 
to opacity. 
none 


idet 
Non LTE calculation of 
switch. 
none 


id3 
tbloop option swtich. 


ieinc 
Flag for x-ray deposition. 
0 


ihop 
Flag to bypass all radiation transport 
calculations. A value of 1 implies a hydro only 
simulation. 


ijz 
im 
Number of values of dr(i) to read in. 


iplot 
Flag for plotting abosorption coefficients. 
none 


iprint 
Flag for simulations prints. 


irez 
Flag controlling rezone. A value of zero 
prohibits rezoning. 


irite 
Controls frequency of binary data dumps. 


istop 
Flag to stop simulation after x-ray deposition. 
0 


itmax 
Number of entries in temperature table, 
currently hardwired to 63. 
63 


NO2 


10−6s 


0.5∗ fcng1 


NO2 
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Table III. More RADFLO Input Variables. 


Variable 
Meaning 
Default Value 


jdeb 
Number of bomb cells. 


jradcn 
Radiation vicous heating switch. 


jsmoo 
Controls smoothing of first pulse. 
none 


j1more 
Flag for the nth event. 
none 


korder 
Order of the spline approximation when 
radius-time-fit smoothing is used. 


lmax 
Maximum number of cells +1 


nnax 
Flag that controls photon energy group 
arrays. 
none 


n 
Number of photons needed to strip a 
nitrogen atoms. 


nclmax 
Maximum number of computational 
cycles the simulation is allowed to run. 
10,000 


ndthmx 
Maximum number of hydro cycles per 
radiation time step. 
10 


nknots 
Number of knots used in determining the 
spline fit when radius-time-fit smoothing 
is used. 


2 


rho 
Ambient air density in gm/cc. 
none 


rhofx(jdeb) 
Initial bomb mass densities. 
none 


scathv 
The photon energy above which the 
aborption coefficient is constant, to 
compensate in an approximate way for 
scattering. 


25 keV 


t 
Blackbody temperature (keV) of x-ray 
source. 
none 


tbloop () 
Output times. 
none 


tcut 
Minimum temperature (eV) for which the 
implicit diffusion treatment will be used 
without a check being made of the actual 
optical depth. 


10 


tfcng1 
The time)s) at which fcng1 may be 
changed. 
0.001s 


tfrac 
This constant, multipled by the expected 
minimum time, gives the time at which 
the debris boundary is frozen in 
Langrangian space. 


th 
The second temperature of a two- 
temperature black body source, if desired. 
none 


thiplt 
The earliest time on the power time plots. 


time 
Problem start time, equal to zero for a 
new start. 


timfit 
timlmt 
Simulation end time. 
none 


uc(jdeb) 
Initial bomb velocity profile. 
none 


ucut 
Minimum optical depth for which the 
implicit diffusion treatment will be used if 
the temperature is less than tcut. 


3.0 


yoke 
X-ray energy (kilotons) to be deposited. 
none 


10−5s 
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A nuclear explosion in air at near-sea level densities produces a distinctive fireball 
optical power output vs time signature consisting of two pulses with a deep minimum in 
between.2 (The same statement applies for infrared and near-UV wavelengths.) For explosions 
at a given altitude the time of the power minimum is uniquely related to the bomb yield, and 
the same is true for the second maximum. Thus by measuring the minimum time or the 
second maximum time one can measure the yield. For altitudes above sea level the depth of 
the minimum decreases with increasing altitude until, at about 40 km (but depending on 
yield), the minimum disappears altogether. 
The RADFLO and HYCHEM codes compute, among other things, the 
fireball radiated power (in several wavelength bands) as a function of time, and also the 
fireball visual radius vs time. Table IV summarizes the sea level air density cases that we 
have run—with and without the nonequilibrium chemistry. The quantities tabulated include 
the input device yield and mass and the timings of seven different phenomenological 
features, most of which are evident in the computed optical power vs time curves. Table V 
summarizes the high-altitude cases we have run. In a few cases there are entries both from a 
HYCHEM run and from an equivalent RADFLOs run. 


Figure 7 is a log-log plot of the computed first maximum, minimum, and second 
max times vs the bomb yield for the sea level cases.3 For the HYCHEM runs, the minimum 
and second max time points are fitted very well by a pair of log-log straight lines, except for 
very small yields below 1 kt, where the tmin points fall below the line. The two straight lines 
are given by the equations 


tmin(Si) = 0.00430Y 0.414 seconds 


and 


t2 max(Si) = 0.0450Y 0.441 seconds. 


2 Note that the optical power (in Watts) is approximately equal to the product of the fireball surface area times 
the optical brightness (in Watts/cm2 steradian) timesπ. So the power and the brightness (or luminosity) are not 
to be confused. 


3 Here the spectral distribution is chosen as that of an unfiltered silicon photodetector. 
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Yield 
(kt) 
Mass 
(g) 
altitude 
(km) 
(s) 
(s) 
(s) 
(s) 


15. 
5.0E+5 
21. 
- 
6.3E-3 
5.0E-2 
7.E-4 


100. 
4.6E+5 
50. 
- 
- 
5.E-4 
- 


t1max 
tmin 
t2max 
tdbshkc 


Table IV. Summary of computed optical power results, for silicon spectral response, 
for explosions at sea level altitude. This table was constructed from RADFLOz 
simulations using the original hydro algorithm. Twelve debris zones were used for 
each simulation. 


Table V. Summary of computed optical power results, for silicon spectral response, 
for high-altitude explosions. 


Yield 
(kt) 
Mass 
(g) 
(s) 
(s) 
(s) 
(s) 
(s) 
(s) 
(s) 


1. Runs at sea level air density without chemistry. 


0.02 
2.7E+4 
3.E-5 
5.3E-4 
8.5E-3 
- 
2.E-5 
- 
- 


0.10 
4.6E+4 
4.E-5 
1.0E-3 
1.6E-2 
- 
2.5E-5 
- 
- 


0.20 
5.8E+4 
5.E-5 
1.4E-3 
2.2E-2 
- 
3.2E-5 
- 
1.1E-3 


0.30 
6.7E+4 
5.E-5 
2.5E-3 
2.6E-2 
- 
1.9E-5 
2.3E-5 
1.1E-3 


1.0 
1.0E+5 
5.5E-5 
4.5E-3 
4.5E-2 
- 
2.3E-5 
4.4E-5 
1.1E-3 


2.0 
1.3E+5 
7.E-5 
5.5E-3 
6.2E-2 
- 
2.7E-5 
7.0E-5 
1.2E-3 


" 
2.5E+5 
8.E-5 
" 
" 
- 
3.3E-5 
- 
1.4E-3 


" 
5.0E+5 
1.0E-4 
" 
" 
- 
7.E-5 
- 
2.1E-3 


" 
1.0E+6 
1.1E-4 
" 
" 
- 
- 
- 
- 


" 
2.0E+6 
1.2E-4 
" 
" 
- 
- 
- 
- 


" 
4.0E+6 
1.8E-4 
" 
" 
- 
- 
- 
- 


" 
8.0E+6 
2.8E-4 
" 
" 
- 
- 
- 
- 


20.0 
2.7E+5 
2.5E-4 
1.4E-2 
1.7E-1 
- 
4.5E-5 
2.E-4 
2.E-3 


" 
5.0E+5 
2.0E-4 
1.4E-2 
1.7E-1 
- 
4.5E-5 
2.E-4 
2.E-3 


100.0 
4.6E+5 
6.0E-4 
2.9E-2 
3.5E-1 
9.E-6 
5.5E-5 
6.E-4 
9.E-3 


1000.0 
1.0E+6 
1.7E-3 
7.5E-2 
9.5E-1 
3.E-4 
9.5E-5 
1.5E-3 
2.0E-2 


" 
2.0E+6 
1.5E-4 
- 


" 
4.0E+6 
1.1E-3 
7.5E-2 
9.5E-1 
3.E-5 
1.2E-4 
1.E-3 
2.E-2 


5000.0 
5.0E+6 
3.4E-3 
1.5E-1 
6.E-4 
1.5E-4 
2.E-3 
3.E-2 


2. Runs at sea level air density with chemistry. 


20.0 
2.7E+5 
2.5E-4 
1.4E-2 
1.7E-1 
- 
4.E-5 
2.5E-4 
1.8E-3 


100.0 
4.6E+5 
5.5E-4 
2.9E-2 
3.5E-1 
9.E-6 
5.5E-5 
6.E-4 
9.E-3 


1000.0 
1.0E+6 
1.5E-3 
7.5E-2 
1.0E+0 
3.5E-4 
9.5E-5 
1.5E-3 
1.8E-2 


t1max 
tmin 
t2max 
tveil 
t50km/s 
tdbshkc 
t0.8eVbp 
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Yield (kt) 


Max. and Min. Times (s) 



10-1 


100 


10-2 


10-3 


10-4 


10-5 
10-2 
10-1 
100 
101 
102 
103 
104 


RADFLOS: First Maximum, Minimum, Second Maximum 
HYCHEM: First Maximum, Minimum, Second Maximum 


Figure 7. Plots of the times of first maximum, minimum, and second maxi- 
mum (in seconds) as functions of the yield (in kilotons), as computed for 
the spectral response function of an unfiltered silicon photodetector for the 
cases in Table IV. The straight lines are scaling expressions fitted to the 
HYCHEM points (see text). 
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The computed minimum and second max times are not affected by nonequilibrium chemistry 
inclusion or by the bomb mass—at least within the range of bomb mass values that we have 
tried. 
For the cases from RADFLOs both the tmin and t2 max points fall slightly below the 
HYCHEM points, consistent with the observation, noted previously, that the RADFLOs 
opacities in the relevant temperature and quantum energy range are slightly smaller than the 
HYCHEM opacities. Presumably if more RADFLOs cases were run the tmin and t2 max points 
would also fall along a pair of straight lines, lying slightly below the two scaling law lines in 
Fig. 7. 
Although the minimum and second maximum times scale neatly with the yield and 
are independent of the bomb mass, that is not true for the first maximum, as can be seen in 
Fig. 7 and Table IV. And, for the range bomb masses that we have run, seems to increase 
with increasing bomb mass for a yield of 2 kt but decreases with increasing bomb mass 
for yields of 20 kt and 1 Mt. 
Most of the energy radiated by a sea level fireball comes out in the second pulse, and 


the total is roughly 30% of the bomb yield. Since the second maximum time varies as Y0.4, 


it follows that the radiant power at second maximum varies roughly as Y0.6. For large yields 
the computed power is somewhat reduced when chemistry is included. 
As noted earlier, for a given yield the time to first maximum depends on the device 
mass and other design details. Moreover, in these idealized computations, which assume a 
spherically symmetric bomb, the first maximum (in the case of the HYCHEM runs) often 
shows distinct fine-structure features: 


1. On the rise before first maximum, there is usually a shoulder or a subsidiary maximum 
that occurs when the shock velocity, which is falling rapidly, passes through 50 km/ 
s (t50km/s in Table IV). The brightness of a shock is a maximum at that velocity (in sea 
level-density air) (see Zinn and Anderson 1973 and Appendix A). After this point the 
power usually flattens out or decreases slightly. 
2. Soon thereafter the inner debris-driven shock wave catches up with the outer fireball 
shock, causing its velocity and brightness to increase (in Table IV). In many cases this 
time is also the time of maximum power—i.e., the main first maximum. At this catch-up 
the bomb-debris kinetic energy adds to the energy of the outer shock, in effect increasing 
the yield driving the shock. After the first maximum the power falls rapidly until the 
temperature of the shock drops to about 0.8 eV. 
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3. As the shock temperature continues to drop below 0.8 eV, the effective specific heat of 
the air decreases (equation of state of air), giving rise to a momentary increase in the 
shock temperature and an additional small maximum (or inflection) in the radiated power. 


Also, the time of the power rise at the beginning of the first pulse can be strongly 
influenced by the so-called x-ray veil. The veil is a diffusely heated region outside the early- 
time fireball which is produced during the initial bomb x-ray emission by the higher energy 
x-rays whose mean free paths are larger than the fireball. If veil temperatures outside the 
fireball exceed about 0.7 eV, the veil hides the early evolution of the fireball. There is an 
abrupt increase in radiated power when the fireball “breaks through” the x-ray veil (tveil in 
Table IV). The x-ray veil temperatures are a strong function of the device mass and other 
design details. The veil is generally not significant in small-yield or large-mass cases. 
Since real nuclear bombs are rarely spherical, the detailed features in the first thermal 
pulse which depend on device mass, etc., can be expected to be smeared out in an actual 
fireball radiant power vs time signature, in contrast to the highly structured first pulse which 
we compute with the spherical RADFLO code. (Note that the present code is able to compute 
details in the first pulse which earlier versions of the code could not.) To investigate this 
point we ran a number of 2 kt calculations using different values of bomb mass, as enumerated 
in Table IV. It can be noted that the computed first maximum time increases with increasing 
mass, and the fine structure disappears. An actual power-time curve might be expected to 
resemble an average of many RADFLO computed power-time curves with different device 
masses. 
The reason for the departure of the tmin times from the scaling law for small yields 
below about 0.3 kt seems to be related to the fact that the shock temperature at tmin increases 
with decreasing values of the yield, while the temperature at the time of the post-first-max 
decreasing-specific-heat blip decreases with decreasing yield, leading to a coincidence of 
the minimum time and the blip when the yield is about 0.2 kt. When the yield is 0.02 kt the 
blip occurs after the minimum. 
The main fireball radiative and hydrodynamic growth is not affected by 
nonequilibrium chemistry. Thus the timings of the various features in the power-time curves 
are not affected either. However, the radiated power magnitudes can be affected strongly— 
especially at early times, and especially for large-yield cases. In the standard 1 Mt case the 
chemistry reduces the height of the first maximum by a factor of five. 
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As the shot altitude is increased, the time duration of the optical pulse decreases, as 
does the depth of the minimum. It is beyond the scope of this report to discuss bhangmeter 
scaling laws for high-altitude events, but they are discussed in Zinn, Kodis, and Sutherland 
(1974). Another high-altitude phenomenon (at ca 20 km) seems to be a distinct splitting of 
the first maximum pulse into a velfb = 50km / s pulse and a debris-shock catch up pulse. 
A few sample computed power-time plots are shown in Figures 8, 9, and 10. 
Figure 8 shows a pair of power-time plots (done with HYCHEM and RADFLOs, 
respectively) for a 1 kt bomb with a mass of 100 kgm in sea level-density air. By our above- 
described bomb-energy-partition model, the debris-surface radiating temperature came out 
to be only 1.2 eV, and essentially all of the bomb yield ended up as internal and kinetic 
energy (assumed 50-50), with less than .001 kt radiated away (i.e. no x-rays). The initial 
surface expansion velocity was about 260 km/s, and the debris-shock velocity was only 
slightly larger. At that velocity the debris shock drove a strong radiative precursor wave 
ahead of it, which moved faster than the shock itself, and whose effective radiating 
temperature was only about 2 eV—much less than the actual Rankine-Hugoniot temperature 
of the shock (Zinn and Sutherland 1981, 1982). As the shock expands outward its velocity 
decreases rapidly. By 15 microseconds it has dropped to about 70 km/s, at which point the 
separated precursor is beginning to form its own shock, and its brightness increases as the 
velocity decreases. Thus the fireball radiant power increases very rapidly, due to the combined 
effects of the increasing brightness and the increasing fireball surface area. At about 23 
microseconds the power passes through a small maximum, when the shock velocity is 50 
km/s (the velocity at which the shock brightness is maximum). After that the power drops 
momentarily before continuing to rise. At about 44 microseconds the debris shock catches 
up with the separate shock from the precursor, causing a momentary surge in shock strength 
and a corresponding increase in brightness. After this time the shock velocity and the 
luminosity decrease rapidly. The first maximum (i.e. the maximum of luminosity × surface 
area) comes slightly after the shock catch up time. After first maximum both the power and 
the brightness decrease smoothly. (The small discontinuous drop at 106 microseconds is a 
numerical artifact arising from a logical switch to stop using the SSW patch at that point.) 
At about 1 millisecond the shock temperature has dropped to about 0.8 eV—the point at 
which the air-specific heat starts to fall rapidly. This results in a momentary temperature 
increase and a brief surge in power. After that time the shock strength continues to drop, as 
does its luminosity. But at the same time, the shock begins to become transparent to visible 
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t (s) 


Slicon Power (W) 


1013 


1012 


1011 


1010 
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10-4 
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10-2 
10-1 
100 


Figure 8. Computed radiant power (silicon band) vs time for a 1 kiloton 
explosion of a device with mass 100 kgm in sea level-density air. HYCHEM 
(solid line) and RADFLOs (dashed curve) results are shown. Both simulations 
used eight equally massed bomb cells with identical partitioning of x-ray, 
kinetic, and internal energies. 
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light (due to a drop in fractional ionization), and one begins to see the shock and the hotter 
air farther in as separate entities (“shock breakaway”). In this time regime there is a developing 
tradeoff between the decreasing brightness of the shock and its increasing transparency to 
the hotter and brighter air within, so that the radiant power, which has been decreasing, 
begins to increase again. Then we have passed through tmin (about 4.3 milliseconds). 
During this same time period the temperatures in the fireball interior have been 
continually equilibrating by radiation transport (at far-UV wavelengths), leading to a structure 
approximating an “isothermal sphere.” After tmin the fireball thermal (and optical) radiation 
comes mainly from this isothermal sphere, which continues to expand hydrodynamically 
and to cool. The rapidly increasing surface area causes the radiated power to increase for a 
while, but eventually, at t2 max, the effects of the combined radiative and hydrodynamic cooling 
become dominant and the radiated power starts to drop. In this time regime, as the fireball 
temperature drops so does its fractional ionization, so that the fireball becomes increasingly 
transparent and its emissivity drops. Eventually, after about 3 × t2 max, it becomes possible to 
see the more luminous debris-air cloud in the fireball interior. The debris is evident in the 
computed brightness profiles but does not have a noticeable effect on the power-time curves. 
In the RADFLOs computation the minimum and second maximum times are slightly earlier 
than in the HYCHEM run, while the first maximum time is considerably later. The first 
maximum power is also quite a bit lower. 


Figure 9 shows a pair of bhangmeter power-time curves for a 1 Mt bomb with a 
mass of 1000 kgm in sea level-density air. Some of the details are quite different from the 1 
kt results. In this case our bomb energy partition model leads to a debris radiating temperature 
greater than 2 kilovolts, and 86.6% of the yield is radiated away as x-rays, leaving only 67 
kt each in kinetic and internal energy. The debris-surface expansion velocity is about 1500 
km/s. The large x-ray yield and high-radiating temperature lead to a very pronounced x-ray 
veil. The visual radius of the veil is about 56 meters and its luminosity is fairly low while it 
completely obscures the smaller and brighter fireball within. The early fireball, produced 
during the x-ray deposition, is a fairly isothermal sphere with a temperature of about 100 eV 
and a radius of about 40 meters. For the next 60 microseconds or so, the fireball expands 
rapidly by radiation transport, but its expansion velocity drops rapidly as the fireball engulfs 
more air and cools. By 60 microseconds the velocity has dropped to a point where 
hydrodynamic motions can compete and a shock forms at the front. The subsequent expansion 
proceeds mostly by hydrodynamics. However, none of these details are visible from the 


40 



t (s) 


Silicon Power (W) 


1015 


1014 


1013 


1012 


10-5 
10-4 
10-3 
10-2 
10-1 
100 
101 


Figure 9. Computed radiant power (silicon band) vs time for a 1 Megaton 
explosion of a device wiht mass 1000 kgm in sea level-density air. 
HYCHEM (solid curve) and RADFLOs (dashed curve) results are shown. 
Again, both simulations used eight bomb cells with identical partitioning of 
x-ray, kinetic, and internal energies. 
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outside because they occur inside the x-ray veil. The fireball breaks through the veil at 
about 300 microseconds, giving rise to a sharp increase in radiated power. The debris shock 
does not catch up with the outer fireball shock until about 1.5 milliseconds, by which time 
it is quite weak and produces only a minimal effect on the brightness. The first maximum 
occurs shortly after the debris-shock catch up. From this point on the fireball expansion 
phenomenology is qualitatively similar to the 1 kt case. The RADFLOs and HYCHEM 
curves are quite similar in this 1 Mt case, although again the RADFLOs minimum and 
second maximum times are slightly earlier than the corresponding HYCHEM times. The 
veil breakthrough time is slightly later in the RADFLOs run. 
We have also run 1 Mt cases with larger bomb masses. With a mass of 4000 kgm 
the x-ray temperature is lower than in the 1000 kgm case and the x-ray yield fraction is 
reduced to 60.7%, with the result that the x-ray veil is much less pronounced. The veil 
breakthrough time is reduced to 30 microseconds. 
We have also run a 1 Mt sea level case with a mass of 1000 kgm with inclusion of 
chemistry. In this case the times of all features in the power-time signal were the same as 
those in Fig. 9, but the height of the first maximum was reduced by a factor of five, and the 
height of the second maximum was reduced by a factor of two. The height of the minimum 
was unchanged. 
To illustrate the grossly different structure of bhangmeter signals for high-altitude 
events, we show in Figure 10 a computed power-time curve for a 100 kt shot at 50 km 
altitude. The figure shows both HYCHEM and RADFLOs results. In this case there is no 
minimum. At this altitude hydrodynamics plays a minor role compared to radiation transport. 
The fireball becomes transparent and the debris cloud becomes visible at 
0.015 s. The debris radiation does not materially affect the radiant power vs time curves.The 
final radiated yield at 0.1 s is 34% of the original 100 kt in the HYCHEM simulation. 
Another useful fireball observable is its radius vs time evolution. Before the Trinity 
test in early 1945 it was predicted by Bethe et al. (1947) that the early fireball growth 
should be approximated by a similarity solution, known currently as the Taylor-Sedov blast 
wave. The blast-wave model is based on the assumptions that gamma (the heat-capacity 
ratio) is a constant and that no appreciable energy is lost by radiation, and it leads to the 
expectation that the fireball radius should increase as the 2/5 power of the time (see e.g. 
Bethe 1964. and Zel’dovich and Raizer 1967). It is interesting to compare it with RADFLO- 
computed radius vs time results. 
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Figure 10. Computed radiant power in the silicon band vs time for a 100 
kiloton explosion of a device with mass 460 kgm in air at 50 km altitude. The 
HYCHEM (solid line) and RADFLOS (dashed line) each used eight debris 
cells, Gudunov hydrodynamics, and identical initial energy partitioning. 
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A plot of the computed fireball radius vs time for 1 kt, 20 kt, and 1 Mt cases are 
shown in Figure 11. The curves are by no means straight lines on the log-log scales. However 
consider the 1kt case, between the time of debris-shock catch up (when for the first time 
the full yield is communicated to the outer shock) and the time of the 0.8 eV blip (after 
which gamma is increasing rapidly), there is a linear segment that corresponds to a power 


law Rfb = const × t0.36—not too different from the 2/5 power. The 20 kt and 1 Mt cases 


have similar slopes for the linear part of the log-log plot. Note that the quantity plotted 
is the visual fireball radius, defined as the locus of optical depth = 1. Hence the collapse to 
zero as the fireball becomes transparent. 
We believe that the current HYCHEM model represents a significant improvement 
over earlier RADFLO versions, but some persistent weak points are discussed in 
Appendix B. 


6.0 SUMMARY 
This report has described the RADFLO family of computer codes for the modeling 
of nuclear fireballs, including their general characteristics and a historical discussion of 
how they have evolved since the original RADFLO code of 1963. We have described the 
similarities and differences between the current LANL codes RADFLOz and HYCHEM 
and the updated MRC version RADFLOs. We have made extensive improvements in the 
codes during the past few years, including (1) replacement of the original Richtmyer-von 
Neumann hydrodynamics algorithm with a one-dimensional version of the modern Godunov 
Riemann-solver algorithm, (2) extension of our steady-state shock and thermal wave (SSW) 
computations to include lower ambient air densities for high-altitude applications, and (3) 


improvements in the method of applying the SSW patch correction for computing fireball 
luminosity and radiant power at early times. We have described these new developments 
in some detail, including a brief discussion of non-equilibrium chemistry effects and their 
computation in HYCHEM. 
We have also made a large number of test computations with the improved codes. 
The results are described, including scaling laws for the times of first maximum, minimum, 
and second maximum as functions of the yield, and variations due to bomb mass, shot 
altitude, and departures from chemical equilibrium. For a few typical cases we have included 
a qualitative description of the phenomenology and growth of the fireball and observables 
such as radiant power vs time and radius vs time. Some comparisons between RADFLOz 
and RADFLOs results are also included. 
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Figure 11. Computed fireball radius vs time for a 1 kiloton, 20 kiloton, and 1 
Megaton explosions, from botton to top respectively. The default parameters 
were used for producing bomb masses of 100, 270, and 1000 kgm, respec- 
tively, in eight debris cells. The simulations were done using HYCHEM and 
Godunov hydro. This is the visual fireball radius, defined as the locus of 
optical depth = 1. Therefore the fireball radius collapses to zero as the 
fireball becomes transparent. 
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Appendix A: SSW Computations 
To extend the range of HYCHEM to altitudes above sea level a set of SSW 
computations were run covering an extended range of ambient air densities. The range 
now extends from downward to . These computations were based on the same algorithms 
and same computer code as described in Zinn and Sutherland (1981). The results are 
summarized in Figure 12, which shows computed values of effective radiating temperature 
for visible light plotted vs the front propagation velocity, for six different values of the 
background air density. 


Appendix B: Numerical Artifacts 
The use of a finite-difference grid requires that material velocities and state variables 
are averaged over grid cells as an approximation to physical reality. It is hoped that the 
approximation improves as cell dimensions are reduced. We have done a large number of 
numerical experiments to examine this point, but are always limited by available computer 
time, since the running time is proportional to the −2 power of the cell thickness 
(approximately). At all levels of cell fineness we see some details in the computed fireball 
front velocity that we do not fully understand, and these details lead to 
corresponding variations in fireball brightness and power. The velocity is computed by 
numerical differentiation of the fireball radius vs time—an inherently noisy process. A 
bothersome question is how well or badly do we compute the process of transition between 
an expanding thermal-wave front and the shock wave that develops from it. That process 
evolves on a spatial scale much smaller than our mesh size, so that we can not resolve it. 
However it occurs when the front velocity is in the range 70 to 100 km/s, when the brightness 
is relatively low but increasing and the power is on the steep rise toward first maximum 


(unless the fireball is still obscured by the x-ray veil). Therefore the details of the radiative 
expansion to shock-wave transition do not affect the power-time curves very much. 
We have already described our procedure for patching the separately computed SSW 
brightness results into the radiation transport computation when the front velocity exceeds a 
threshold value vquitfit , which is usually set at 30 km/s. We discontinue implementation 
of the patch when the velocity drops below that value, and use the standard RADFLO 
radiation transport algorithm to compute the brightness. In some cases this change in the 
algorithm results in a small discontinuous drop or jump up in the computed power. Figure 


8 shows such a discontinuity at 1.06 ×10−4s. 
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Figure 12. Plots of effective radiating temperature (for visible light) vs the 
wavefront propagation velocity for six different values of the background air 
density. (Open circle–density=0.00123 g/cm3; solid circle–density=0.0001; solid 
square–density=0.0005; open square–density=0.0002; down triangle–den- 
sity=0.0001; up triangle–density=5×10-5.) 
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An area of much-needed improvement in the HYCHEM model is the treatment of 
gamma rays and neutrons. They are approximated in the current model by a pair of analytic 
expressions, which may be wrong in some cases by an order of magnitude. Given today’s 
easily accessible computational power, deposition via a real radiation transport algorithm 
seems feasible. We believe this would vastly improve the pre-first-maximum part of the 
power-time curve. Second, the bomb debris could be modelled as a second turbulent fluid 
because of the debris/air interface, Taylor unstable nature. The influence of bomb debris on 
the free electron density and the total opacity could then be directly calculated. 
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